clear 
est clear

********************************************************************************
*	Programs
******************************************************************************** 

*** program for the OLS regressions... 
capture program drop regs
program define regs, eclass

	*** Syntax 
	syntax varlist(numeric) [, treatment(varlist numeric) rounds(numlist) ctrls(varlist numeric)]
 
	*** loop over outcomes and rounds
	qui foreach y of local varlist {
		foreach r in `rounds' {
		    
			*** create the block variables, always centred for sample with nonmissing outcome 
			capture drop blocks*
			qui tab block, gen(blocks)
			qui foreach v of varlist blocks* {
				replace `v' = . 			if missing(`y') & round ==`r'
				su `v' 						if round == `r' & !missing(`y')
				replace `v' = `v' - r(mean) if round == `r' & !missing(`y')
			}
			
			*** Model, includes: 
			*	- treatment,
			*	- block FE fully interacted,
			*	- clustered SE's at club level,
			*	- tracking weights (=1 in RR and 12m fups)
			* 	- controls fully interacted (when adjusted)
			local model_name = "`y'_`r'"
			* if (length("`ctrls'")>0) local model_name = "`y'_`r'_adj"
			
			eststo `model_name': reg `y' b0.`treatment'##c.(`ctrls' blocks*) [pw=w] if round==`r', vce(cluster club)
			
			*** save stats 
			* control mean and SD 
			qui su `y' [aw=w] if `treatment'==0 & round==`r' 
			estadd scalar c_mean = r(mean) 	
			estadd scalar c_sd = r(sd) 
			* test coeffs p-values, when treatments are combined
			if ("`treatment'"=="treat_ther") { 
				* p-value for IPT-G
				test 1.treat_ther == 0 
				estadd scalar pval_ther_any = r(p)
			}
			* p-value for the two treatments and comparison stats 
			else if ("`treatment'"=="treat") {
			    * p-value for IPT-G
				test 1.treat == 0 
				estadd scalar pval_ther = r(p)
				* p-value for IPT-G+
				test 2.treat == 0 
				estadd scalar pval_cash = r(p) 
				* diff between treaments 
				lincom 2.treat - 1.treat
				estadd scalar est_diff = r(estimate) 
				estadd scalar pval_diff = r(p) 
			}
		}
	}
	drop blocks* 
end


*** program to extract q values for the outcomes and treatment specified, within each round
capture program drop get_qvals
program define get_qvals, eclass

	*** Syntax 
	syntax varlist(numeric) [, rounds(numlist) tname_in(string) tsuff_out(string) method(string)]
	
	*** Create vector of p-values
	* matrix to store them
	local varcount : word count `varlist' 
	matrix pvals_tab = J(`varcount',1,.)  
	matrix colnames pvals_tab = pvals
	
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`r'
			est replay
			* collect pvals in matrix 
			if ("`tname_in'" == "est_diff") {
				lincom 2.treat - 1.treat 
				matrix pvals_tab[`i',1] = `r(p)'
				matrix list pvals_tab				
			}
			else {
				matrix list r(table)
				matrix pvals_tab[`i',1] = r(table)["pvalue","`tname_in'"]
				matrix list pvals_tab			    
			}
		}
		* save in a dataset 
		capture drop pvals`r'
		svmat pvals_tab, names(col)
		rename pvals pvals`r'
		
		*** get q-values  
		capture drop qvals`r'
		qqvalue pvals`r', method(`method') qvalue(qvals`r') 
	}
	
	*** add to estimates 
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`r' 
			* add as a scalar
			local qv = qvals`r' in `i'
			di `qv'
			estadd scalar qval_`tsuff_out' = `qv'
		}
	} 
end




********************************************************************************
*	I - Unadjusted results in appendix
********************************************************************************

use "${data}/Analysis_wide.dta", clear
distinct block club


*** prep the data 
* get it in long format
reshape long phq8_score ghq12_score rosb_score res_score loc_score ghq12_min phq8_min w, i(block cr_id club treat ) j(round)
* keep only what we need 
keep block cr_id club treat treat_ther round w ///
	 age_b everpregnant_b evermarried_b PPI_score_b phq8_score_b ///
	 phq8_score ghq12_score rosb_score res_score loc_score ghq12_min phq8_min
* drop missings, use the phq score that was always asked of everyone
drop if missing(phq8_score)


*** Regressions
* RR: use the combined IPT-G treatment indicator 
regs phq8_min phq8_score ghq12_min ghq12_score rosb_score res_score loc_score, ///
	treatment(treat_ther) rounds(1) 
* 12m: skip the missing outcomes
regs phq8_min phq8_score ghq12_min ghq12_score , ///
	treatment(treat) rounds(2) 
* 24m: everything
regs phq8_min phq8_score ghq12_min ghq12_score rosb_score res_score loc_score, ///
	treatment(treat) rounds(3) 

	
*** MHT adjustments 
* primary outcomes 
get_qvals phq8_min ghq12_min , rounds(1) tname_in("1.treat_ther") tsuff_out("ther_any") method("simes")
get_qvals phq8_min ghq12_min , rounds(2 3) tname_in("1.treat") tsuff_out("ther") method("simes")
get_qvals phq8_min ghq12_min , rounds(2 3) tname_in("2.treat") tsuff_out("cash") method("simes")
* secondary outcomes 
get_qvals phq8_score ghq12_score rosb_score res_score loc_score, rounds(1) tname_in("1.treat_ther") tsuff_out("ther_any") method("simes")
get_qvals phq8_score ghq12_score , rounds(2) tname_in("1.treat") tsuff_out("ther") method("simes")
get_qvals phq8_score ghq12_score , rounds(2) tname_in("2.treat") tsuff_out("cash") method("simes")
get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(3) tname_in("1.treat") tsuff_out("ther") method("simes")
get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(3) tname_in("2.treat") tsuff_out("cash") method("simes")
* difference between the two coefs
get_qvals phq8_min ghq12_min , rounds(2 3) tname_in("est_diff") tsuff_out("diff")  method("simes")
get_qvals phq8_score ghq12_score , rounds(2) tname_in("est_diff") tsuff_out("diff") method("simes")
get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(3) tname_in("est_diff") tsuff_out("diff") method("simes")


*** Unadjusted table 
* RR:
#delimit ;
esttab 	phq8_min_1 ghq12_min_1 phq8_score_1 ghq12_score_1 rosb_score_1 res_score_1 loc_score_1
	using "${tables}/results-mh-revised-1.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval_ther_any qval_ther_any , fmt(3 3 0 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* 12m:
#delimit ;
esttab 	phq8_min_2 ghq12_min_2 phq8_score_2 ghq12_score_2 
	using "${tables}/results-mh-revised-2.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat 2.treat) 
	stats(c_mean c_sd N pval_ther qval_ther pval_cash qval_cash est_diff pval_diff qval_diff, fmt(3 3 0 3 3 3 3) 
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values" "Coeff.: IPT-G+ - IPT-G" "H0: IPT-G=IPT-G+ p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* 24m: 
#delimit ;
esttab 	phq8_min_3 ghq12_min_3 phq8_score_3 ghq12_score_3 rosb_score_3 res_score_3 loc_score_3
	using "${tables}/results-mh-revised-3.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat 2.treat) 
	stats(c_mean c_sd N pval_ther qval_ther pval_cash qval_cash est_diff pval_diff qval_diff, fmt(3 3 0 3 3 3 3) 
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values" "Coeff.: IPT-G+ - IPT-G" "H0: IPT-G=IPT-G+ p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 







est clear
********************************************************************************
*	II - Adjusted results for appendix
********************************************************************************

use "${data}/Analysis_wide.dta", clear
distinct block club


*** prep the data 
* get it in long format
reshape long phq8_score ghq12_score rosb_score res_score loc_score ghq12_min phq8_min w, i(block cr_id club treat ) j(round)
* keep only what we need 
keep block cr_id club treat treat_ther round w ///
	 age_b everpregnant_b evermarried_b PPI_score_b phq8_score_b ///
	 phq8_score ghq12_score rosb_score res_score loc_score ghq12_min phq8_min
* drop missings, use the phq score that was always asked of everyone
drop if missing(phq8_score)


*** Regressions
* RR: use the combined IPT-G treatment indicator 
regs phq8_min phq8_score ghq12_min ghq12_score rosb_score res_score loc_score, ///
	treatment(treat_ther) rounds(1) ctrls(phq8_score_b age_b everpregnant_b evermarried_b PPI_score_b)
* 12m: skip the missing outcomes
regs phq8_min phq8_score ghq12_min ghq12_score , ///
	treatment(treat) rounds(2) ctrls(phq8_score_b age_b everpregnant_b evermarried_b PPI_score_b)
* 24m: everything
regs phq8_min phq8_score ghq12_min ghq12_score rosb_score res_score loc_score, ///
	treatment(treat) rounds(3) ctrls(phq8_score_b age_b everpregnant_b evermarried_b PPI_score_b)

	
*** MHT adjustments 
* primary outcomes 
get_qvals phq8_min ghq12_min , rounds(1) tname_in("1.treat_ther") tsuff_out("ther_any") method("simes")
get_qvals phq8_min ghq12_min , rounds(2 3) tname_in("1.treat") tsuff_out("ther") method("simes")
get_qvals phq8_min ghq12_min , rounds(2 3) tname_in("2.treat") tsuff_out("cash") method("simes")
* secondary outcomes 
get_qvals phq8_score ghq12_score rosb_score res_score loc_score, rounds(1) tname_in("1.treat_ther") tsuff_out("ther_any") method("simes")
get_qvals phq8_score ghq12_score , rounds(2) tname_in("1.treat") tsuff_out("ther") method("simes")
get_qvals phq8_score ghq12_score , rounds(2) tname_in("2.treat") tsuff_out("cash") method("simes")
get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(3) tname_in("1.treat") tsuff_out("ther") method("simes")
get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(3) tname_in("2.treat") tsuff_out("cash") method("simes")
* difference between the two coefs
get_qvals phq8_min ghq12_min , rounds(2 3) tname_in("est_diff") tsuff_out("diff")  method("simes")
get_qvals phq8_score ghq12_score , rounds(2) tname_in("est_diff") tsuff_out("diff") method("simes")
get_qvals phq8_score ghq12_score rosb_score res_score loc_score , rounds(3) tname_in("est_diff") tsuff_out("diff") method("simes")


*** Unadjusted table 
* RR:
#delimit ;
esttab 	phq8_min_1 ghq12_min_1 phq8_score_1 ghq12_score_1 rosb_score_1 res_score_1 loc_score_1
	using "${tables}/results-mh-revised-1-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N pval_ther_any qval_ther_any , fmt(3 3 0 3)  
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* 12m:
#delimit ;
esttab 	phq8_min_2 ghq12_min_2 phq8_score_2 ghq12_score_2 
	using "${tables}/results-mh-revised-2-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat 2.treat) 
	stats(c_mean c_sd N pval_ther qval_ther pval_cash qval_cash est_diff pval_diff qval_diff, fmt(3 3 0 3 3 3 3) 
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values" "Coeff.: IPT-G+ - IPT-G" "H0: IPT-G=IPT-G+ p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* 24m: 
#delimit ;
esttab 	phq8_min_3 ghq12_min_3 phq8_score_3 ghq12_score_3 rosb_score_3 res_score_3 loc_score_3
	using "${tables}/results-mh-revised-3-adj.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat 2.treat) 
	stats(c_mean c_sd N pval_ther qval_ther pval_cash qval_cash est_diff pval_diff qval_diff, fmt(3 3 0 3 3 3 3) 
	label("Control mean" "Control SD" "Observations" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values" "Coeff.: IPT-G+ - IPT-G" "H0: IPT-G=IPT-G+ p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 